%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute an off-diagonal block of the matrix discretizing a single contour.
% Note that this INCLUDES corrections near the diagonal for handling the
% singularity in the kernel.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = get_A_single_offd(C,ind1,ind2,kh)

ntot     = size(C,2);

[Y_g1,  X_g1  ] = meshgrid(C(1,ind2), C(1,ind1));
[Y_g2,  X_g2  ] = meshgrid(C(4,ind2), C(4,ind1));
[Y_dg1, X_dg1 ] = meshgrid(C(2,ind2), C(2,ind1));
[Y_dg2, X_dg2 ] = meshgrid(C(5,ind2), C(5,ind1));

ima = sqrt(-1);
nn1 = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
nn2 = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
dd  = sqrt((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2);

K   = funcK(kh,nn1,nn2,(Y_g1 - X_g1),(Y_g2 - X_g2),dd);

A   = K.*sqrt(Y_dg1.^2 + Y_dg2.^2);

% Correction for the singularity
MU6 = [ 0.2051970990601250e1 + 0.2915391987686505e1;...
        -0.7407035584542865e1 - 0.8797979464048396e1;...
        0.1219590847580216e2 + 0.1365562914252423e2;...
        -0.1064623987147282e2 - 0.1157975479644601e2;...
        0.4799117710681772e1 + 0.5130987287355766e1;...
        -0.8837770983721025   - 0.9342187797694916];
D = abs(ind1'*ones(1,length(ind2)) - ones(length(ind1),1)*ind2);
D = abs(mod(D + 6, ntot) - 6);
[ii1, ii2] = find(D.*(D<7));
for j = 1:length(ii1)
    i1 = ii1(j);
    i2 = ii2(j);
    d  = D(i1,i2);
    A(i1,i2) = (1 + MU6(d))*A(i1,i2);
end

return
